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The construction and operation of large scale quantum information devices presents a 
grand challenge. A major issue is the effective control of coherent evolution, which requires 
accurate knowledge of the system dynamics that may vary from device to device. We 
review strategies for obtaining such knowledge from minimal initial resources and in an 
efficient manner, and apply these to the problem of characterization of a qubit embedded 
into a larger state manifold, made tractable by exploiting prior structural knowledge. We 
also investigate adaptive sampling for estimation of multiple parameters. 
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1. Introduction 

Recently, much effort has been put into the construction of large scale quantum 
devices operating in the coherent regime. This has been spurred by the possibilities 
offered by quantum communication and information processing, from secure 
transmission, simulation of quantum dynamics, to the solution of currently 
intractable mathematical problems [1]. Many different physical systems have been 
proposed as the basic architecture upon which to construct quantum devices, 
ranging from atoms, ions, photos, quantum dots and superconductors. For large 
scale commercial application, it is likely that this will involve scalable engineered 
and constructed devices with tailored dynamics requiring precision control. 

Due to inevitable manufacturing tolerances and variation, each device will 
display different behaviour even though they may be nominally "identical". For 
their operation, they will need to be characterized as to their basic properties, 
response to control fields, and noise or decoherencc [2]. Wc may also need to know 
how ideal is the system in the first place, for instance the effective Hilbert space; 
what we may assume to be a qubit may have dynamics involving more than two 
effective levels. Extracting this information efficiently and robustly is crucial. 

In a laboratory setting, an experimentalist may have access to many tools with 
which to study a system, e.g. spectroscopy and external probes. In a production 
setting, provision of these extra resources may be difficult, expensive, or impossible 
to integrate with the device. It is therefore an important to understand what sort 
of characterisation can be performed simply using what is available in situ. Ideally, 
we would also like to be able to characterize the performance of a device with as 
little prior knowledge of its behaviour, e.g. how it responds to control fields, if this 
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is the information which we are trying to obtain in the first place. Characterization 
using only the in situ resources of state preparation and measurement, even where 
it is possible, is challenging, due to the increasing complexity of the signals, 
number of parameters to be estimated and the complexity of reconstructing a valid 
Hamiltonian from the resulting signal parameters. Robust and efficient methods 
of data gathering and analysis, preferably as automated as possible, are therefore 
essential. Here, we review progress in tackling these problems and we present new 
results for the problem of characterizing the Hamiltonian of a qubit embedded in 
a larger state manifold. 



A standard tool for determining discrete quantum dynamics is quantum process 
tomography (QPT). Assume a completely positive trace preserving map A acting 
on a system state p by 



where Xj are the Kraus operators satisfying XjXj = 1. QPT is a method by 
which we can determine {A^} by seeing how initial states of the system evolve 
under the map. For a complete characterisation, a complete set of initial states 
{pj} and quantum state tomography on the final states are both required ^. 
In general for a d-dimensional system, we need to use cF — 1 input states, and 
quantum tomography either requires a measurement with at least cP outcomes, for 
instance, symmetric informationally complete positive operator- valued measures 
SIC-POVMs [5], or projective measurements in several different bases. 

The ability to generate a complete set of initial states is a strong assumption. 
In many physical systems, it is only possible to directly prepare a set of orthogonal 
states, e.g. by projective measurement, or even only a single "fiducial" state. The 
usual assumption that we can generate any state by coherent control of a fiducial 
state is invalid when we are trying to determine how to control the system in 
the first place. Equally, most systems can only be measured directly in a single 
fixed basis, and other measurement bases are assumed to be available through 
coherent control. This shows that QPT cannot be used as a starting for quantum 
system identification in a setting where control has not been already established. 
We require a mechanism by which we can bootstrap our knowledge and abilities 
until control can be enacted upon the system. 

Assuming the system dynamics can be approximated to a reasonable degree 
by Hamiltonian dynamics, the first core challenge is the identification of the 
intrinsic system, Hamiltonian H = Ej\Ej){Ej\, which can be specified by the 
energy eigenstates and eigenvalues (up to rescaling). Optimum protocols for 
identification of the Hamiltonian dynamics depend on the available resources. One 
general paradigm introduced in [6] assumes a situation where we are restricted to 
preparation and measurement in a single fixed basis {|ej)}, i.e., we can prepare 

^ There are variants of QPT which use ancillas, entanghng or other coherent operations which 
may offer certain advantages [3, 4]. However they all require operational resources beyond what 
we assume for system characterisation. 
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the system in one of the basis states \ej), allow it to evolve for some time t 
before projectively measuring in the same basis, obtaining state \ek) with some 
probability Pjk{t). The computational basis in this case can be defined in terms 
of preparation and measurement and our task is to characterize the system 
Hamiltonian in this basis, i.e. obtain the eigenenergies Ej and the eigenstates 
\Ej) (up to a global phase), solely from the data traces of Pjk{t). 

Hamiltonian characterization in this paradigm has been considered in a number 
of papers. It has been shown that a generic Hamiltonian for a single qubit can 
be recovered from preparation and measurement in a fixed basis up to a certain 
set of phases and an unobservable global energy shift. The extra phases become 
relevant only when additional resources are available that allow us to initialize the 
system in non-measurement basis states, or apply control that alters the system 
Hamiltonian. In the latter setting is was also shown how the relative phases for two 
Hamiltonians could be recovered by composite rotations, vaguely reminiscent of 
Ramsey spectroscopy [7] . If the Hamiltonians do not commute or coincide with the 
preparation/measurement basis, full control over a single qubit can be realized and 
the system and control Hamiltonians can be characterized up to a sign factor [8] . 
The results can be extended to multi-qubit or higher level systems [9]. 

Upon characterisation of the intrinsic system Hamiltonian the effect of 
applying controls must be identified. We can model this by assuming a 
Hamiltonian H{\) that depends on classical control field parameters A. In the 
simplest case we may approximate the response of the system by H{X) = Hq + 

^jHj, where Hq is the intrinsic system Hamiltonian and Hj are perturbations 
resulting from the application of control Xj. This has been considered in [6, 8] 
where the effect of multiple control fields on a single qubit were characterized 
with respect to a reference Hamiltonian. For coherent operation, incoherent effects 
should be small but they will still need to be characterized. Although complete 
characterization of the dynamics for open systems is a daunting task under 
certain simplifying assumptions on the type of decoherence, e.g., pure dephasing 
or relaxation in a natural basis, the number of parameters to be estimated can be 
reduced and the relevant information extracted [10, 11, 12, 13]. 

Although the general characterization paradigm described is quite restrictive, 
in some cases further restrictions must be imposed to deal with limited resources. 
Characterization of the coupling constants in a spin network when only a subset 
of the spins at the boundary can be measured and initialized is an example [14, 
15, 16]. Another is estimation of leakage out of a subspace or coupling to unknown 
states, where we generally cannot measure the populations of these states directly. 
General bounds on subspace leakage were derived in [17] based on the Fourier 
spectrum of the observed Rabi oscillations. This approach is useful when there 
are potentially many states very weakly coupled to the subspace of interest. In 
many cases, however, leakage may be due to coupling to a small number of states 
outside the subspace, e.g., when we encode a qubit using the lowest two states 
of a slightly anharmonic oscillator. In this case a control applied to the qubit 
transition will induce some coupling to the third (nuisance) level, giving rise to 
unwanted dynamics. Characterization of this coupling allows the design of pulses 
that can suppress such unwanted excitations. 
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3. Hamiltonian estimation for embedded qubit 

Formally, we consider a three-level system subject to a control field resonant with 
the 1—2 transition (see also [18]), where |1) and |2) can be regarded as the qubit 
states. Assuming a constant amplitude field f{t) = Acos{ujt), transforming to a 
rotating frame and making the rotating wave approximation, this leads to an 
effective Hamiltonian of the form 




H=\di d2 . (3.1) 



We choose a rotating frame where the off-diagonal elements dn are real and 
positive. If higher order processes such as two-photon transitions between states 1 
and 3 are negligible, we can assume ds = 0. The objective is to characterize both 
the qubit transition coupling di as well as the coupling to the nuisance level c?2 
and the detuning 5 as a result of the anharmonicity. 

If the system can be initialized in the basis states \n) for ri = 1,2, 3 and 
we perform complete projective measurements in the same basis, then the 
probabilities 

can be determined for all k,i for different times tj. However, if we can only 
initialize and measure the system in the ground state then only a single population 
evolution trace pii{t) is available. For a two-level system this is sufficient to infer 
the population of p22 but the presence of the third level means that pu{t) only 
gives us limited information about the population of the other levels. 

The existence of a non-zero detuning 5 complicates the problem substantially. 
In the absence of anharmonicity, i.e., for 6 = analytic expressions for pii{t) 
can be obtained and the problem reduced to a single frequency estimation 
problem [19]. For 5^0 there is no simple closed form for the signal pii{t) 
and the eigenvalues are no longer of the form 0, ±A, but are instead LO12 = 
A2 — Ai = a; — Aoj and LO23 = A3 — A2 = w -|- Aco. For small detunings 6 relative 
to the coupling strengths dn, the frequency splitting Aco is much smaller than 
Lo. To obtain obtain a frequency resolution of Au using spectral analysis would 
require a signal length of at least l/(Aa;). However, using the structure of the 
signal and restricting to solutions consistent with our prior knowledge we can do 
considerably better. We know that the observed signal must be of the form 

Pii{t) = ao + ai cos((a; — Auj)t) + a2 cos((a; -|- Auj)t) -f 03 cos(2a;t). (3.3) 

A natural starting point for a maximum likelihood estimation is thus to 
choose basis functions go = 1, gi{t) = cos((a; — Auj)t), g2{t) = cos((a; -|- Auj)t) and 
gsit) = cos(2c(jt), and following standard techniques, maximize the log-likelihood 
function [9, 20] 



L{{uj, Au}\d) oc - logio 



mb(h^)'' 



(3.4) 
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Figure 1: "Qubit" with leakage to a slightly detuned third level. The raw signal 
pu{t) (low-descrepancy sampling, A't = 100, A'^e = 100) provides (top, left) a 
preliminary estimate of the main frequency. The DFT (bottom left) shows only 
a single peak. The log-likelihood (right) clearly shows that a model with a small 
splitting of Aw around 0.1 (white ellipse) is significantly more likely than a single 
frequency model. (Colour online) 

where = 4 is the number of basis functions, Nt is the number of data points 
and 

n=0 m=0 

The elements hm of the {rrib, l)-vector h are projections of the (1, Nt)-data vector 
d onto a set of orthonormal basis vectors derived from the non-orthogonal basis 
functions gm{t) evaluated at the respective sample times t„- Concretely, setting 
Gmn = gm{tn), let Am and Bm be the eigenvalues and corresponding (normalized) 
eigenvectors of the x mi, matrix GG"^ with G= {Gmn), and let E= (e^'m) 
be a matrix whose columns are e^i- Then we have H = VG and h = i?dt with 
y = diag(am^'^^)£^^^, and the corresponding coefficient vector is ai = h.^V [9]. 

Fig. 1 shows that the log-likelihood function of the data provides strong 
evidence for a non-zero detuning, even though no peak splitting is detectable 
in the Fourier spectrum of the signal. For the given input data the log-likelihood 
function has a squeezed peak that is narrow in one direction but much broader 
in the other. The fact that the peak is squeezed in a direction not aligned with 
a coordinate axis shows that the uncertainties in oj and Acj are not independent. 
The plot also shows that the width of the peak along the Auj direction is much 
greater than that in u direction. 
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Figure 2: Standard deviation of LOest and Awest for 256 simulated experiments each 
as a function of the number of time samples Nt and the number of experiment 
repetitions N^.. The scale is logj^o- (Colour online) 

This is also reflected in the observed uncertainties of the estimates of a; and 
Acj. If we take the (wgst, Awgst) to be the coordinates for which the log-likelihood 
peaks then their standard deviations give an indication of the uncertainty in our 
estimates. Fig. 2 shows the standard deviation of ujest and Aa;est for 256 simulated 
experiments each as a function of the number of time samples Nt and the number 
of experiment repetitions Nf.^ on a logarithmic scale. The plots look qualitatively 
similar, suggesting a similar scaling, but the scale shows that the uncertainty 
of the Aw estimates is about one order of magnitude greater than that of the 
a; estimates. We observe a similar scaling for the estimated amplitudes Om for 
m = 0, 1, 2, 3 (not shown). The median relative error in the estimated Hamiltonian 
shows a similar behaviour though with some kinks, and it is interesting to note 
that the relative errors in the Hamiltonian tend to be larger than the errors in 
the estimated frequencies and amplitudes of the signal. 

Preliminary results suggest that fragility of the reconstruction procedure is 
responsible for the observed larger spread in the errors of the reconstructed 
Hamiltonian even if the uncertainty of the estimated signal parameters is quite 
low as shown in Fig. 3 (left). The reconstruction procedure can sometimes fail, 
leading to outliers in the relative Hamiltonian error histogram shown in Fig. 3 
(right), which typically correspond to unphysical Hamiltonians. 

If the likelihood function does not have a sufficiently sharp peak, we can 
adaptively refine the sampling to reduce the uncertainty. A simple adaptive 
strategy is as follows: 

1. Preliminary sampling: Over a pre-chosen sampling period, measure at 
randomly chosen sampling times and compute the likelihood of various 
models based on this initial data. 

2. Uncertainty estimate: The uncertainty of the solution can be estimated from 
the sharpness and relative height of the highest peak in the likelihood plot. 
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Figure 3: (left) Error scaling as a function of sampling times and accuracy, (right) 
Histogram of reconstruction error showing outliers corresponding to unphysical 
solutions. (Colour online) 

3. Refinement: From the initial data, choose an ensemble of probable models 
and calculate the weighted expected variance of their data traces as a 
function of time. 

4. New samples: We make additional measurements at times for which the most 
probable models differ the most and use the new data points to update the 
estimates of the signal parameters and Hamiltonian. 

5. Repeat as necessary. 

Although similar in spirit to the adaptive Bayesian identification strategy 
proposed Wiseman et al. for single parameter Hamiltonian estimation [21], we 
do not minimize the variance of a single parameter as the Hamiltonian depends 
on multiple parameters. Another difficulty is that the multi-parameter likelihood 
function is usually far from Gaussian and the expectation values of uj and Aw 
tend to differ substantially from the maximum likelihood estimate. In this case 
using the expectation values and variances with respect to a given parameter is 
not necessarily a good indicator of the real uncertainty of the model. 

We apply this to the above qutrit system. Starting with a 100-point low- 
discrepancy sampling of the selected time range [0,20], we obtain A'^g = 100 
measurements per time point. The resulting log-likelihood function, shown in 
Fig. 4(left), has a squeezed peak centred at (w, Aw) = (1.9468, 0.1087). To narrow 
the peak width we resample using the initial estimate for (a;,Aa;). A simple 
and computationally cheap strategy is to choose new sample points at integer 
multiples of T/2 = 'K/ujf>st-, and get a few accurate samples using e.g. = 1000. 
The idea is that we will be most sensitive to small modulations in the peak 
heights due to Aw at these times. Indeed we find that the relative errors for the 
frequency and amplitude estimates improve substantially. Fig. 4 (right). Yet, the 
relative error of the reconstructed Hamiltonians does not always follow the same 
trend and sometimes actually increases. Further analysis shows that this is due 
to the reconstruction step and the complex dependence of H on the estimated 
parameters. 
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log Likelihood(u,ALj) — original data log Likelihood(u,Aij) — with new data 




Figure 4: Log-likelihood comparision (left) Initial, (right) After adaptive sampling. 
The probability is more highly concentrated and the accuracy of the solution 
greatly improved. (Colour online) 



This suggests that it would be desirable to estimate the Hamiltonian 
paramaters directly from the data, maximizing the likelihood 

-Nt/2 



P(d|{r2, a, e}) oc exp 



Nt 



(3.6) 



or its logarithm, and for convenience we have defined di = Oleosa, ^2 =r2sina, 
and 8 = 4:6. Although there is no simple closed form for pii{tj) in this case, it 
can be computed numerically and the challenge is finding the global optimum 
of the 3-parameter likelihood function. Without any prior information, we 
first compute the log-likelihood on a 3D grid of parameter values, find the 
region where the global optimum is expected, and use this information as 
a starting point of a local optimization routine to find the peak. For the 
example above, this yields (fi, a, e)opt = (1-7159, 0.9494, 0.5340) for the initial 
data, versus (fi, a, e)opt = (1.7323, 0.9557, 0.4999) with the new data (actual 
(1.7321,0.9553,0.5000)). The relative error of the reconstructed Hamiltonian 
for the inital data is 0.0363, comparable (even slightly higher) to the estimate 
obtained using the previous two-step approach. Using the new data, the relative 
error is substantially lower, 3.8672 x 10~^ if we maximize (3.6) instead of 
maximizing (3.4) followed by reconstruction. This approach appears suitable for 
systematic adaptive estimation, which will be investigated further in future work. 



4. Conclusion 

We have outlined the problem of system characterisation and why it is an 
essential basic building block of quantum control for many prospective quantum 
information processing devices. We have shown that such a bootstrapping 
procedure is possible and that much information can be gained through utilisation 
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of limited in situ initially present operational capabilities. By exploiting prior 
knowledge and reasonable assumptions on the structure and behaviour of the 
system, maximum likelihood analysis offers efficient and robust characterisation 
and reconstruction of complex systems. 

Scalability remains a challenge. In the general setting, slightly increasing the 
size of the system leads to an explosion in signal complexity [9] and directly 
applying system identification to three of more qubits is a formidable task. 
Some complexity reduction can be achieved using Bayesian signal estimation 
in order to split up frequency and amplitude estimation, though this has some 
drawbacks in ensuring physically allowed reconstructed Hamiltonians. Extending 
the Bayesian estimation directly to Hamiltonians would alleviate reconstruction 
validity but direct optimisation of the likelihood is challenging for more than a 
few parameters. Hence there is the need for a similar complexity reduction for 
Hamiltonian parameter estimation. Exploiting as much structural information 
about the system is essential in making the problem tractable, especially in the 
case of restricted resources. 

Adaptive estimation for multiple parameters is a ripe area for further 
exploration. Practical online schemes may require pragmatic methods of 
determining adaptive measurements, as full Bayesian optimisation involves 
integration over a highly peaked payoff function in a high dimensional parameter 
space. Development of effective yet computationally efficient optimisation routines 
is imperative. 

Relaxing more assumptions or expanding the set of resources available would 
give greater experimental relevance. Preparation and measurement capabilities 
may vary, for instance instead of projective measurement, some experimental 
proposals implement continuous [22], weak or generalized measurements. States 
may also be prepared by relaxation or adiabatic passage and may not coincide 
with the measurement basis. 

More generally, there are interesting connections between the compressive 
sensing (CS) and sparse reconstruction paradigm [23, 24] and how our model- 
based systen characterisation techniques work. Instead of the union of low- 
dimensional (linear) sub-spaces model in CS, we instead have a solution space 
as the union of low-dimensional manifolds of parameters. An extension of the 
notion of "basis incoherence" and general techniques for efficient reconstruction 
from sparse data would be extremely beneficial. 
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